#####################################################
#
# Python2 script to reproduce Fig. 1 in the comment
# authors: Thomas Lang, Stephan Hesselmann
#
#####################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import glob

ecolor = ['#3e77af','#60af92','#efc600','#cc2909','#3e77af','#60af92','#efc600','#cc2909']
fcolor = ['#ffffff','#ffffff','#ffffff','#ffffff','#a3c1e5','#99d1b9','#eddea2','#ea6868']
ecolor = ['#cc2909','#efc600','#60af92','#3e77af','#cc2909','#efc600','#60af92','#3e77af']
fcolor = ['#ffffff','#ffffff','#ffffff','#ffffff','#ea6868','#eddea2','#99d1b9','#a3c1e5']
marker = ['o','s','D','^','o','s','D','^']
markerSize = np.array([5,4.5,4,5.5,5,4.5,4,5.5])*2.0

fsTicksLabel = 18
fsTicksLabelInset = 16
fsAxesLabel = 21
fsLegend = 16
fsLabel = 18

a1 = np.array([np.sqrt(3.), 0.])
a2 = np.array([np.sqrt(3.)/2., 3./2.])
b1 = np.array([2.*np.pi/np.sqrt(3.), -2.*np.pi/3.])
b2 = np.array([0., 4.*np.pi/3.])
KK = (2*b2 + b1)/3.
K = { 6 : np.array([4.1887902, 0.]), 9 : np.array([4.1887902, 0.]), 12 : np.array([4.1887902, 0.]), 15 : np.array([4.1887902, 0.]) }
def E(k):
    ek = np.zeros(len(k))
    for ik in range(len(k)):
        ek[ik] = np.abs(1. + np.exp(1.j * np.dot(k[ik], -a2)) + np.exp(1.j * np.dot(k[ik], a1 - a2)))
    return ek


i_L = 0
i_gamma = 1
i_U = 2
i_q = 3
i_kx = 4
i_ky = 5
i_delta = 6
i_sigma = 7

L_list = [ 6, 9, 12, 15 ]
gamma_list = [ 0. ]
U_list = [ 0.5, 3.75 ]

fig = plt.figure(figsize=(11*0.7,8*0.7), dpi=96)
ax = fig.add_subplot(111)
ax.tick_params(which='both', axis='both', direction='in', bottom=True, top=True, left=True, right=True, labelsize=fsTicksLabel)
ax.set_xlabel(r"$a \Delta k$", fontsize=fsAxesLabel)
ax.set_ylabel(r"$E / t$", fontsize=fsAxesLabel)

ax.set_xlim(-0.07, 2.05)
ax.set_ylim(-0.07, )
plt.axhline(0, color='#dddddd', zorder=-2, linewidth=0.5)
plt.axvline(0, color='#dddddd', zorder=-2, linewidth=0.5)
ax.xaxis.set_minor_locator(ticker.AutoMinorLocator(2))
ax.set_yticks([0, 0.4, 0.8, 1.2, 1.6])
ax.yaxis.set_minor_locator(ticker.AutoMinorLocator(2))
ax.text(0.51, 0.62, r"$v_0$", color='#888888', transform=ax.transAxes, fontsize=fsLabel, ha='left')
ax.text(0.51, 0.31, r"40% reduced $v_0$", color='#cc2909', transform=ax.transAxes, fontsize=fsLabel, ha='left')
ax.text(0.85, 0.10, r'$\gamma = 0$', color='#000000', transform=ax.transAxes, fontSize=fsLabel, ha='left', 
			bbox=dict(boxstyle="round",
                   ec=(219/255., 202/255., 175/255.),
                   fc=(244/255., 232/255., 211/255.),
                   )
         )

filelist = glob.glob("data.txt")
filelist.sort()

k = np.zeros((400,2))
kabs = np.zeros(400)
for ik in range(len(k)):
    k[ik,:] = KK + np.sqrt(3) * KK[1] * (np.pi*(b2)/np.linalg.norm(b2)) * ik/len(k)
    kabs[ik] = np.linalg.norm( KK[1] * (np.pi*(b2)/np.linalg.norm(b2)) * ik/len(k))*3.

ax.plot(kabs, E(k)*1., color='#aaaaaa', ls=":", linewidth=1.2)
k = np.linspace(0., 1.14, 2)
ax.plot(k, 3.**0.5/2.*0.6 * k, color='#cc2909', ls="-", linewidth=1.2)
ax.plot(k, 3.**0.5/2. * k, color='#aaaaaa', ls="-", linewidth=1.4)

cnt_U = 0
cnt_L = 0
for filename in filelist:
    with open(filename) as f:
        lines = (line for line in f if not line.startswith('L'))
        new_del = []
        for l in lines:
            new_del.append(','.join(l.split()))
        data = np.loadtxt(new_del, delimiter=',', skiprows=0)

    data_filter = list(filter(lambda x : x[i_L] in L_list and x[i_gamma] in gamma_list and x[i_U] in U_list, data))

    ip = 0
    p = {}
    for U in U_list:
        data_filter_U = np.array(list(filter(lambda x : x[i_U] == U, data_filter)))
        for L in L_list:
            data_filter_L = np.array(list(filter(lambda x : x[i_L] == L, data_filter_U)))

            q_list = data_filter_L[:, i_q]
            kx_list = data_filter_L[:, i_kx]
            ky_list = data_filter_L[:, i_ky]
            delta_list = data_filter_L[:, i_delta]
            sigma_list = data_filter_L[:, i_sigma]
                        
            if U == 0.5:
                label = r"$U/t="+str(U)+",\, \ L="+str(L)+"$"
            else:
                label = r"$U/t="+str(U)+", L="+str(L)+"$"
            if U > 0.5:
            	sc = 0.6
            else:
            	sc = 1.0
            p[ip] = ax.errorbar(q_list, delta_list, yerr=sigma_list, markeredgewidth=1.6, capsize=8, fmt=marker[cnt_L], markersize=sc*markerSize[cnt_L], color=ecolor[cnt_L], ecolor=ecolor[cnt_L], mfc=fcolor[cnt_L], label=label, clip_on=False)
            ip += 1
            if L == 15 and U == 0.5:
                ax.plot(q_list[:2], delta_list[:2], '--', color='#cc2909', linewidth=1.4)
            if L == 15 and U == 3.75:
                ax.plot(q_list[:2], delta_list[:2], '-', color='#cc2909', linewidth=1.4)
            cnt_L += 1
        cnt_U += 1

# The published version of this plot includes the linear finite size extrapolations
# provided by Tang et al.; The data for was provided in pdf form only, 
# which we merged with the pdf generated from this script.

plt.legend(loc=(0.7, 0.01), prop={'size': fsLegend}, frameon=False, labelspacing=0.4, handletextpad=-0.25)
l1 = plt.legend([p[0], p[1], p[2], p[3]], [r'$U/t=0.5,\, \ L=6$', r'$U/t=0.5,\, \ L=9$', r'$U/t=0.5,\, \ L=12$', r'$U/t=0.5,\, \ L=15$'], loc=(0.02, 0.7), prop={'size': fsLegend}, frameon=False, labelspacing=0.35, handletextpad=0.25, borderpad=0.25, facecolor='wheat', framealpha=0.5)
l2 = plt.legend([p[4], p[5], p[6], p[7]], [r'$U/t=3.75, L=6$', r'$U/t=3.75, L=9$', r'$U/t=3.75, L=12$', r'$U/t=3.75, L=15$'], loc=(0.02, 0.42), prop={'size': fsLegend}, frameon=False, labelspacing=0.35, handletextpad=0.25, borderpad=0.25, facecolor='wheat', framealpha=0.5)
plt.gca().add_artist(l1)    

plt.tight_layout()
plt.savefig('Fig1.pdf')
